N = 32
h = 1 / N;
disp('Iterative method using FFT');
tic
[G, V] = speciter2d(N);
toc
disp('Direct method using HSC');
tic
Tree = five_setup(V, N, h);
f = eye(N*N);
Gfive = five_solve(N,Tree,f);
toc
disp('relative 2-norm difference');
norm(G-Gfive,2)./norm(G,2)

Ind = reshape(1:N*N,N,N);
SourceY = 1:N/4;
SourceX = N/2+1:3*N/4;
IndSource = reshape(Ind(SourceY,SourceX),[],1);

TargetY = 1:N/4;
TargetX = 1:N/4;
IndTarget = reshape(Ind(TargetY,TargetX),[],1);

Sspec = svd(G(IndTarget, IndSource));
Sfive = svd(Gfive(IndTarget, IndSource));
rspec = length(find(Sspec/max(Sspec)>1e-4));
rfive   = length(find(Sfive  /max(Sspec)>1e-4));
disp('spectral method rank')
rspec
disp('finite difference rank')
rfive

